home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
TeX 1995 July
/
TeX CD-ROM July 1995 (Disc 1)(Walnut Creek)(1995).ISO
/
web
/
fweb
/
fweb-1.40
/
demos
/
adj.web
< prev
next >
Wrap
Text File
|
1993-10-29
|
66KB
|
2,278 lines
% $Id: adj.web,v 4.9 90/06/16 07:32:04 karney Exp $
% $Log: adj.web,v $
% Revision 4.9 90/06/16 07:32:04 karney
% More documentation.
%
% Revision 4.8 90/06/15 13:59:56 karney
% Inserted more user documentation on how to compile and run the code.
%
% Revision 4.7 90/05/11 17:37:26 karney
% Cleaned up use of macros to provied generic UNIX version. This should
% probably be restructure to provide more orthogonality.
%
% Revision 4.6 90/05/07 17:42:25 karney
% Fix typo in readline addition.
%
% Revision 4.5 90/05/07 12:24:15 karney
% Use readline library to get terminal input on Sun
%
% Revision 4.4 90/04/26 16:22:13 karney
% Second attempt to get version number
%
% Revision 4.3 90/04/26 16:19:54 karney
% Include version number in title
%
% Revision 4.2 90/04/19 14:46:35 karney
% Minor cosmetic changes following John's suggestions
%
% Revision 4.1 90/04/19 12:11:43 karney
% Initial entry of Aug 22, 1989 version of code
\def\title{Adjoint function for a toroidal plasma; $ $Revision: 4.9 $ $}
% Change some of FWEB's formatting macros.
\def\\#1{\leavevmode\hbox{\def\${{\sl \char`\$}}\def\%{\char`\%}%
\it#1\/\kern.05em}} % Italic type for identifiers
\def\@@#1{\leavevmode\hbox{\def\${\char`\$}\def\%{\char`\%}%
\rm#1\kern.05em}} % Roman intrinsic functions
\let\Csc=\rm % Font for comments
% Some general math macros
\def\abs#1{\left\vertbar#1\right\vertbar} % Absolute values
\def\frac#1#2{{#1\over#2}}
\let\d=\partial
\def\pd#1#2{\mathchoice{\d#1\over\d#2}{\d#1/\d#2}{\d#1/\d#2}{\d#1/\d#2}}
% Some shortcuts for this program
\def\u{{\bf u}}\def\S{{\bf S}}
\def\subb#1#2#3{_{#1[#2]#3}}
\def\ja#1{\tilde\jmath_{#1}}
\def\jn#1#2#3{\bar\jmath\subb{#1}{#2}{#3}}
\def\qa#1{\tilde q_{#1}}
\def\s{*}
\def\lave#1{\overline{#1}}
@n
@* Introduction. This program is supposed to calculate the Green's
function for rf-driven current in an axisymmetric tokamak. This is the
solution to the Spitzer-H\"arm problem with the loop voltage $\Phi$ set to
$TL/Q$. The rf-driven current density is given by
$$J_{\parallel 0}= \int d^3\u_0\, \lambda \S_0\cdot\nabla_0\chi,$$
where $\S_0$ is the bounce-averaged rf-induced flux (at the midplane).
In this program, velocities are normalized to an arbitrary reference
velocity $u_n$ (which is usually either $\sqrt{T/m}$ or $c$), times to
$\tau_n=(\Gamma^{e/e}/u_n^3)^{-1}$, current densities to $qnu_n$, fluxes to
$n/(\tau_n u_n^2)$, and $\chi$ to $qu_n\tau_n$.
(1) May 3, 1989:
Initial version with only the first Legendre harmonic running for RF
conference, May 1--3, 1989.
(2) May 29, 1989:
Include an arbitrary number Legendre harmonics.
(3) May 31, 1989:
Convert to using dynamic memory allocation.
(4) August 22, 1989:
Allow choice on boundary conditions at $u_{\rm max}$.
@ This program is written to run on VAX/VMS systems, Sun Sparcstations
(using double precision) and on Crays with cft77 (using single precision).
This is controlled by four macros |VAX|, |SUN|, |CRAY|, |UNIX|. It should be
relatively straightforward to modify to program to run in other
environments. The major system dependent aspect of this code is the use of
dynamic memory. However this is easily turned off using the macros
|STATIC|. Then all memory is taken from a single statically allocated
array of size |MEMSIZE|.
For the VAX, the prgram should be created with:
{\parskip=0pt\parindent=10em\tt
ftangle -mVAX adj\par
fortran/g\_float adj\par
link adj\par}
For the Cray, use:
{\parskip=0pt\parindent=10em\tt
ftangle -mCRAY -d adj\par
{\rm (send {\tt adj.f} to Cray)}\par
cf77 -o adj adj.f\par}
(The {\tt -d} converts |do|/|end do| constructions into labeled |do|
loops.)
For the Sun, use:
{\parskip=0pt\parindent=10em\tt
ftangle -mSUN -=c=\#.web.c adj\par
f77 -I/usr/gnu/include adj.f adj.web.c -o adj $\backslash$\par
\hskip3em -L/usr/gnu/lib -lreadline -ltermcap\par}
If you need to use static memory, then specify {\tt -mSTATIC -m"MEMSIZE
100000"} to ftangle; this uses static memory of size 100000.
@ Input parameters. When you run this program it asks you to input the
following parameters: |t|, |c2|, |ze|, |zi|, |eps|, |nmax|, |umax|, |imax|,
|kmax|, |dt|, |tmax|, |alpha|, |rho|, |lmax|.
|t| is the temperature in units of $m u_n^2$. Thus if you want velocities
normalized to $\sqrt{T/m}$, set $t = 1$. Alternatively if normalization is
to speed of light, set $t = T/(mc^2)$.
|c2| is the square of the speed of light (normalized to $u_n^2$). If you
want to run nonrelativistic cases, then set $t=1$ and ${\it c2} = 10^{20}$.
|ze| is the electron charge state. This is usually 1, but you can get the
Lorentz limit $Z_i\to\infty$ with $|ze|=0$ and $|zi|=1$.
|zi| is the ion charge state.
|eps| is the inverse aspect ratio $r/R$. This is only used in the magnetic
field routine |deltab|. This routine is currently set up to define a
circular cross-section tokamak. $\epsilon=0$ corresponds to a uniform
magnetic field (no trapped particles).
|nmax| is the number of grid points in the $u_0$ direction.
|umax| is the maximum value of $\abs{{\bf u}_0}$. The $u_0$ grid is evenly
spaced between $0$ and $u_{\rm max}$.
|imax| is the number of grid points in the $\theta_0$ direction. The
$\theta_0$ grid is evenly space between $0$ and $\theta_{\rm tr}$, the
trapping angle. $\chi({\bf u}_0)$ is odd about $\theta = \pi/2$ and is
zero for $\theta_{\rm tr} < \theta < \pi - \theta_{\rm tr}$ in the trapped
zone.
|kmax| is the number of integration points for field line integrals. These
are only done once at the beginning of a run, so |kmax| should be set quite
large (e.g., 1000).
|dt| is the time step. This needs to be chosen large enough to get to the
steady state quickly but not so large that the iteration becomes unstable.
A typical value is 1.
|tmax| is the number of time steps, typically 100--300.
|alpha| is a weighting factor in the time stepping algorithm. $\alpha=0$
gives an explicit algorithm. $\alpha=0.5$ gives Crank-Nicholson iteration.
$\alpha=0.55$ (i.e., a small amount of forward weighting) seems to be best.
|rho| governs the treatment of the boundary at $u=u_{\rm max}$. For
$\rho=1$, we have $\chi''(u_{\rm max})=0$, and for $\rho=0$, we have
$\chi'(u_{\rm max})=0$. The boundary has the least effect on the solution
in the interior for $\rho=1$.
|lmax| is the maximum number of Legendre harmonics. $|lmax|=0$ corresponds
to collision off a Maxwellian. For $\epsilon=0$, set $|lmax|=1$ (since
$\chi$ consists of the first Legendre harmonic only). Otherwise try
$|lmax|=21$.
The program cycles through repeatedly asking for input until the end of the
input file (for terminal input this is Control-Z on the VAX or Control-D on
the Sun) or else nonpositive values are read in for |t| or |c2|.
@ Output. This program outputs a small amount of data to the terminal plus
a file of unformatted data.
The terminal output consists of an echo of the input parameters. Every 10
time steps, a line is written out consisting of |time|, |cond|, |abserr|,
|relerr|. |time| is the time step. |cond| is the present estimate of the
electrical conductivity (which this program calculates because $\chi$ is
just the Spitzer-H\"arm function). |abserr| and |relerr| are estimates of
the absolute and relative errors in $\chi$ based on the residue.
The data written to the file consists of the input variables, the $u$ and
$\theta_0$ arrays and a two-dimensional array for $\chi$. See the section
``Write out results to disk.'' Note that $\chi$ and $f_m$ are defined at
the mid-points of the cells i.e., at $u=|ua(n)|$, $\theta_0=|th0(i)|$ for
$0\le n<n_{\rm max}$, $0\le i<i_{\rm max}$.
@ Files. The program talks to three different files, which we call
|input|, |output|, |adjout|. |input| and |output| are normally terminal
input and output. |adjout| is the output file of unformatted data (see
above).
On the VAX, the system supplies a {\tt .dat} file type to the |adjout| file
(i.e., it is normally the {\tt adjout.dat} in the current directory). You
can redirect |input| and |output| from/to a file and rename the file
|adjout| with VAX logical names. For example:
{\parskip=0pt\parindent=10em\tt
define/user sys\$output adj2.in\par
define/user sys\$input adj2.out\par
define/user adjout adj2.dat\par
r adj\par}
On the Sun, you can use normal file redirection for |input| and |output|. You have to rename |adjout| explicity. For example:
{\parskip=0pt\parindent=10em\tt
adj <adj2.in >adj2.out\par
mv adjout adj2.dat\par}
On the Cray, your can use command line arguments to specify file names.
For example:
{\parskip=0pt\parindent=10em\tt
adj input=adj2.in output=adj2.out adjout=adj2.dat\par}
@ References:
[1] B. J. Braams and C. F. F. Karney, {\it Conductivity of a relativistic
plasma}, Phys.~Fluids {\bf 1B}, 1355--1368 (1989).
[2] C. F. F. Karney, N. J. Fisch, and A. H. Reiman,
{\it Green's function for rf-driven current in a toroidal plasma},
in {\it Radio-Frequency Power in Plasmas},
Proc.\ Eighth Topical Conf. on Radio Frequency Plasma Heating,
University of California, Irvine, CA,
May 1--3, 1989,
edited by R. McWilliams,
AIP Conference Proceedings, No. 190
(AIP, New York, 1989), pp.~430--433.
[3] C. F. F. Karney and N. J. Fisch, {\it Efficiency of Current Drive by
Fast Waves}, Phys.~Fluids {\bf 28}, 116--126 (1985).
[4] C. F. F. Karney, {\it Fokker--Planck and Quasilinear Codes},
Computer Physics Reports {\bf 4}, 183--244 (1986).
@ Macros for machine dependencies.
@#if defined(CRAY)
@m CRAY 1
@m SUN 0
@m VAX 0
@m UNIX 0
@#elif defined(SUN)
@m CRAY 0
@m SUN 1
@m VAX 0
@m UNIX 0
@#elif defined(VAX)
@m CRAY 0
@m SUN 0
@m VAX 1
@m UNIX 0
@#else
@m CRAY 0
@m SUN 0
@m VAX 0
@m UNIX 1
@#endif
@ Define precision.
@#if VAX || SUN || UNIX /* Use double precision on the VAX and Sun*/
@m real double precision /* Default precision */
@m areal DBLE /* Coerce to default precision */
@m const(x) x##d0 /* Constant in default precision */
@m epsilon 1.0d-16 /* Approximate round-off error */
@m single REAL /* Coerce to single precision */
@m memunit 8 /* Size of a real in bytes (the VAX's unit of
memory) */
@#endif
@#if CRAY /* Use single precision on the Cray */
@m real REAL
@m areal REAL
@m const(x) x##e0
@m epsilon 1.0e-16
@m single REAL
@m memunit 1 /* Size of a real in words (the Cray's unit of
memory) */
@#endif
@ I/O units. Also set things up to use the readline library if it's
available (to permit editing and recall of input lines).
@m outunit 6
@m diskunit 20
@#if !defined(READLINE)
@#if SUN
@m READLINE 1
@#else
@m READLINE 0
@#endif
@#endif
@#if READLINE
@m LINELEN 80
@m inunit line(1:min(linelen,LINELEN))
@#else
@m inunit 5
@#endif
@ Label definitions.
@m Start #:0
@m Exit #:0
@ The Sun uses |implicit undefined(a-z)| instead of |implicit none|
@f implicit_none implicit
@#if SUN
@m implicit_none implicit undefined(a-z)
@#else
@m implicit_none implicit none
@#endif
@ Abbreviated names for some identifiers. CFT cannot handle identifiers
longer than 8 characters.
@#if CRAY
@m maxwellian maxwel
@m potential pot
@m isotropic isotrop
@m matrixinit matinit
@m legendrej legj
@m legendrejs legjs
@m legendrejn legjn
@#endif
@* The main program.
@a
program adjoint
implicit_none
integer nmax,imax,kmax,tmax,lmax,tskip
real t,c2,ze,zi,eps,umax,dt,alpha,rho
real cpu,cpu0
external operinit,second,runa
@#if READLINE
character line*LINELEN
integer readline,linelen
external readline
@#endif
call operinit /* Open files */
Start: continue
write(outunit,*)
$ 'Enter t,c2,ze,zi,eps,nmax,umax,imax,kmax,dt,tmax,alpha,rho,lmax'
@#if READLINE
linelen=readline("> ",line)
if (linelen<=0) then
write(outunit,*)
goto Exit
end if
@#endif
read(inunit,*,end=Exit) t,c2,ze,zi,eps,nmax,umax,imax,kmax,
$ dt,tmax,alpha,rho,lmax
if (t<=0 || c2<=0) goto Exit
write(outunit,*) t,c2,ze,zi,eps,nmax,umax,imax,kmax,
$ dt,tmax,alpha,rho,lmax
tskip=10
call second(cpu0)
call runa(t,c2,ze,zi,eps,nmax,umax,imax,kmax,dt,tmax,
$ alpha,rho,lmax,tskip)
call second(cpu)
write(outunit,'(1x,a,f10.2,a)') 'CPU time: ',cpu-cpu0,' secs'
goto Start
Exit: continue
stop
end
@<Functions and subroutines@>
@#if SUN || UNIX
@c
@<C functions@>
@#endif
@ Open files for the CRAY.
@m STRING(s) #s /* Macros get unit numbers into argument of |link|. */
@m S(s) STRING(s)
@<Functions...@>=
@#if CRAY
subroutine operinit
implicit_none
external link
call link('input=terminal,unit'//S(inunit)//'=(input,open,text),'// @|
$ 'output=terminal,unit'//S(outunit)//'=(output,create,text),'// @|
$ 'print'//S(outunit)//','// @|
$ 'unit'//S(diskunit)//'=(adjout,create,seq)//')
return
end
@#endif
@ Open files for the VAX. Define a function to get the CPU time
@<Functions...@>=
@#if VAX
subroutine operinit
implicit_none
external lib$init_timer
call lib$init_timer
open(unit=diskunit,file='adjout',form='unformatted',status='new')
return
end
subroutine second(cpu)
implicit_none
integer cputime
real cpu
external lib$stat_timer
call lib$stat_timer(2,cputime) /* |cputime| is in 10 ms increments */
cpu=const(0.01)*cputime /* convert to seconds */
return
end
@#endif
@ Open files for the SUN. Define a function to get the CPU time
@<Functions...@>=
@#if SUN || UNIX
subroutine operinit
implicit_none
integer cputime,clock
external clock
cputime=clock()
open(unit=diskunit,file='adjout',form='unformatted',status='unknown')
return
end
subroutine second(cpu)
implicit_none
integer cputime,clock
real cpu
cputime=clock()
cpu=const(0.000001)*cputime /* convert to seconds */
return
end
@#endif
@ A readline function for the Sun. This reads in a line with editing and
recall using the Gnu readline library.
@<C functions@>=
@#if READLINE
#include <stdio.h>
#include <strings.h>
#include <readline/readline.h>
/* Read a string with a prompt returning string in |result|. String length
is returned as function value. */
long int readline_ (prompt, result, lprompt, lresult)
char prompt[], result[];
long int lprompt, lresult;
{
char *cprompt, *cresult;
cprompt = (char *)malloc (lprompt+1);
strncpy (cprompt,prompt,lprompt);
cprompt[lprompt] = 0;
/* Get a line from the user. */
cresult = readline (cprompt);
free (cprompt);
/* If the line has any text in it, save it on the history. */
if (cresult && *cresult)
add_history (cresult);
if (cresult == (char *)NULL) {
return (EOF);
} else {
strncpy (result, cresult, lresult);
lresult = strlen (cresult);
free (cresult);
return (lresult);
}
}
@#endif
@ Some C compatibility routines for the Sun.
@<C functions@>=
@#if SUN || UNIX
long clock_() /* Interface to C clock routine */
{
return(clock());
}
@#endif
@* Memory allocation. All memory for arrays is allocated in |runa|. When
tangled for static memory allocation, the memory is taken from a single
statically allocated array |memory| of size |MEMSIZE|. For dynamic memory
allocation, memory is obtained using |lib$get_vm| or |mzalloc|. These
return an address of a memory block. This is associated with an array with
a pointer statement (for the Cray) or by passing the address by value down
to |runb|. (See the definition of the macro |mem|.)
@#if !UNIX
@#if defined(STATIC)
@m STATIC 1 /* Allocate memory statically? */
@#else
@m STATIC 0
@#endif
@#else
@m STATIC 1
@#endif
@#if !defined(MEMSIZE)
@m MEMSIZE 4000
@#endif
@ Utility macros for memory allocation. For each array used in |runb|, a
index is defined by prepending a |p| to the name of the array. This index
is then defined using |alloc|. Arrays which are local to a single
subroutine are allocated from a single array |work|. |workalloc| figures
out the size of this array.
@m ptr(x) p##x /* A pointer into |memory| formed by prepending |p| */
@m alloc(x,y) ptr(x)=memsize; memsize=memsize+y;
@m workalloc(y) worksize=max(worksize,y)
@#if STATIC
@m mem(x) memory(ptr(x))
@#else
@#if VAX
@m mem(x) %val(baseaddr+memunit*ptr(x))
@#endif
@#if CRAY || SUN
@m mem(x) memory(ptr(x))
@#endif
@#endif
@ Since we deal only in harmonic numbers $l$ equal to 0, 1, 3, 5, etc., we
also define macros to map between these and an index |la| equation to 0, 1,
2, 3, etc. We also need macros to index into two non-rectangular arrays
|garr| and |harr|.
@m laf(l) (l+1)/2 /* Map |l| onto collapsed index |la| */
@m lf(la) max(2*la-1,0) /* Map |la| onto expanded index |l| */
@m gind(la,la1) ((((la+3)-1)*(la+3))/2 + (la1+1))
@m gv(la,la1) garr(gind(la,la1))
@m hind(la,la1,la2) (((la-1)*la*(2*la-1))/6 + (la1-1)*la + la2)
@m hv(la,la1,la2) harr(hind(la,la1,la2))
@ Allocate the necessary memory and call |runb| to do the actual
computation.
@f @<Declare memory allocation stuff@> integer
@f @<Index declarations@> integer
@f pointer integer
@<Functions...@>=
subroutine runa(t,c2,ze,zi,eps,nmax,umax,imax,kmax,dt,tmax,
$ alpha,rho,lmax,tskip)
implicit_none
integer nmax,imax,kmax,tmax,lmax,tskip
real t,c2,ze,zi,eps,umax,dt,alpha,rho
integer nmax1,imax1,lmaxa,lmaxa1
@<Declare memory allocation stuff@>
external runb
/* |nmax1| and |imax1| are the allocated dimensions for two-dimensional
arrays. These might differ from |nmax| and |imax| on the Crays, where even
strides are inefficiently handled. */
nmax1=nmax
imax1=imax
@#if CRAY /* To avoid memory conflicts on the Crays */
if (mod(nmax1,2)==1) nmax1=nmax1+1 /* Make even for accessing
|ijy(0:nmaxa1,0:6,0:1)| by row in |potential| */
if (mod(imax1,2)==0) imax1=imax1+1 /* Make odd for accessing
|chi(0:imax1-1,0:nmax-1)| by row in |decomp| */
@#endif
/* |lmaxa1| is the allocated dimension for the |la| index. This is
usually |lmaxa| but must be at least $1$ to ensure that dimension
specifications such as |1:lmaxa1| are legal. Also $H_{1,1,1}$ is needed
even when |lmaxa=0|. */
if (mod(lmax,2)==0) lmax=lmax-1 /* Round |lmax| down to an odd number */
lmaxa=laf(lmax)
lmaxa1=lmaxa
if (lmaxa1==0) lmaxa1=1 /* But |lmaxa1| needs to be at least 1 */
memsize=0
worksize=0
@<Index specifications@> /* Define all the array pointers and
figure |worksize| */
alloc(work,worksize)
@<Allocate memory@>
call runb(t,c2,ze,zi,eps,nmax,umax,imax,kmax,dt,tmax,
$ alpha,rho,lmax,tskip,
$ nmax1,imax1,lmaxa,lmaxa1,worksize,
$ mem(ua),mem(ug),mem(th0a),mem(th0g),mem(csth0a),mem(chi),
$ mem(legjy),mem(chil),mem(legp),mem(duug),mem(fu),mem(dtt),
$ mem(fm),mem(psid),mem(lam1a),mem(lam2g),mem(harr),mem(uinv),
$ mem(th0inv),mem(resid),mem(resida),mem(work))
@<Deallocate memory@>
return
end
@
@<Declare memory allocation stuff@>=
integer memsize,worksize,ptr(work)
@<Index declarations@>
@#if STATIC
real memory
dimension memory(0:MEMSIZE-1)
@#else
@#if VAX
integer baseaddr,status,lib$get_vm,lib$free_vm
external lib$get_vm,lib$free_vm
@#endif
@#if CRAY
real memory
dimension memory(0:*)
pointer (baseaddr,memory)
integer status,mzalloc,mzdalloc
external mzalloc,mzdalloc
@#endif
@#if SUN
real memory
dimension memory(0:*)
pointer (baseaddr,memory)
integer status,malloc
external malloc,free
@#endif
@#endif
@ Allocate memory required using the appropriate system calls.
@#if !STATIC
@#if VAX
@m memalloc(size,addr) status=lib$get_vm(size,addr)
@m memdealloc(size,addr) status=lib$free_vm(size,addr)
@m bad(s) mod(s,2)==0
@#endif
@#if CRAY
@m memalloc(size,addr) status=mzalloc(addr,size)
@m memdealloc(size,addr) status=mzdalloc(addr,size)
@m bad(s) s!=0
@#endif
@#if SUN
@m memalloc(size,addr) addr=malloc(size); status=addr
@m memdealloc(size,addr) call free(size)
@m bad(s) s==0
@#endif
@#endif
@<Allocate memory@>=
@#if STATIC
if (memsize>MEMSIZE) then
write(outunit,*) 'Memory needed/allocated: ', memsize, '/', MEMSIZE
return
end if
@#else
memalloc(memunit*memsize,baseaddr)
if (bad(status)) then
write(outunit,*) 'Couldn''t allocate memory: ',memsize
return
end if
@#endif
@
@<Deallocate memory@>=
@#if !STATIC
memdealloc(memunit*memsize,baseaddr)
@#endif
@* Do a case.
@f @<Global declarations@> integer
@<Functions...@>=
subroutine runb(t,c2,ze,zi,eps,nmax,umax,imax,kmax,dt,tmax,
$ alpha,rho,lmax,tskip,
$ nmax1,imax1,lmaxa,lmaxa1,worksize,
$ ua,ug,th0a,th0g,csth0a,chi,legjy,chil,legp,duug,fu,dtt,
$ fm,psid,lam1a,lam2g,harr,uinv,th0inv,resid,resida,work)
implicit_none
integer worksize
real work
dimension work(0:worksize-1)
@<Global declarations@>
@<Initialization@>
do time=0,tmax
@<Make a step@>
end do
@<Write out results@>
return
end
@* The velocity space meshes.
@<Glob...@>=
integer nmax,imax,nmax1,imax1,n,i
real du,umax,ua,ug,dth0,th0max,th0a,th0g,csth0a,chi @/
/* The velocity mesh */
dimension ua(0:nmax-1),ug(0:nmax),th0a(0:imax-1),th0g(0:imax),
$ csth0a(0:imax-1),chi(0:imax1-1,0:nmax-1)
real pi
@ Define the indices for velocity space meshes in main memory of |runa|.
@<Index decl...@>=
integer ptr(ua),ptr(ug),ptr(th0a),ptr(th0g),ptr(csth0a),ptr(chi)
@ Allocate the space needed for these arrays.
@<Index spec...@>=
alloc(ua,nmax)
alloc(ug,nmax+1)
alloc(th0a,imax)
alloc(th0g,imax+1)
alloc(csth0a,imax)
alloc(chi,imax1*nmax)
@ Initialize velocity meshes, $\chi$, etc.
@<Init...@>=
du=umax/nmax
do n=0,nmax-1
ua(n)=(n+const(0.5))*du
ug(n)=n*du
end do
ug(nmax)=umax
pi=4*atan(const(1.0))
th0max=atan2(const(1.0),sqrt(deltab(pi,eps))) /* $\sin^{-1}(1/\sqrt b)$*/
dth0=th0max/imax
do i=0,imax-1
th0a(i)=(i+const(0.5))*dth0
csth0a(i)=cos(th0a(i))
th0g(i)=i*dth0
end do
th0g(imax)=th0max
do n=0,nmax-1
do i=0,imax-1
chi(i,n)=0
end do
end do
@* Potentials. The following routines calculate the potentials and their
derivatives. First define arrays to hold Legendre functions evaluated on
the velocity mesh.
@m jl(n,a,la) legjy(n,a,0,la)
@m yl(n,a,la) legjy(n,a,1,la)
@m djl(n,a,la) legjy(n,a,2,la)
@m dyl(n,a,la) legjy(n,a,3,la)
@<Glob...@>=
real c,c2
real legjy /* Arrays holding Legendre functions */
integer lmax,lmaxa,lmaxa1
dimension legjy(0:nmax1-1,0:6,0:3,0:lmaxa1)
external initjy
@
@<Index decl...@>=
integer ptr(legjy)
@
@<Index spec...@>=
alloc(legjy,nmax1*7*4*(lmaxa1+1))
@ Fill the arrays.
@<Init...@>=
c=sqrt(c2)
lmaxa=laf(lmax)
call initjy(ua,c,nmax,nmax1,lmaxa,lmaxa1,legjy,
$ work(0),work(2*(lf(lmaxa)+2+1)*7),lf(lmaxa)+2)
@ Define |initjy|.
@<Functions...@>=
subroutine initjy(ua,c,nmax,nmax1,lmaxa,lmaxa1,
$ legjy,jarray,djarray,lmax1)
implicit_none
integer nmax,nmax1,lmaxa,lmaxa1,n,mult,la,l,a,lmax1
real ua,legjy,c,jarray,djarray
dimension ua(0:nmax-1),legjy(0:nmax1-1,0:6,0:3,0:lmaxa1)
dimension jarray(-lmax1-1:lmax1,0:6),djarray(-lmax1-1:lmax1,0:6)
external legendrejn
do n=0,nmax-1
call legendrejn(ua(n),c,lf(lmaxa),jarray,djarray,lmax1)
do la=0,lmaxa
l=max(lf(la),0)
mult=(-1)^(l+1)
do a=0,6
jl(n,a,la)=jarray(l,a)
yl(n,a,la)=mult*jarray(-l-1,a)
djl(n,a,la)=djarray(l,a)
dyl(n,a,la)=mult*djarray(-l-1,a)
end do
end do
end do
return
end
@ Figure the workspace for |initjy| (|jarray| and |djarray|).
@<Index spec...@>=
workalloc(2*(lf(lmaxa)+2+1)*7 * 2)
@ Calculate the potentials and their derivates for the $l$th Legendre
harmonic.
These are some definitions so we don't have to pass a whole slew of arrays
around to subroutines.
@m j(n,a) legjy(n,a,0) /* Short-hand expressions for $j$, $y$, their
derivatives and integrals. */
@m y(n,a) legjy(n,a,1)
@m dj(n,a) legjy(n,a,2)
@m dy(n,a) legjy(n,a,3)
@m ij(n,a) ijy(n,a,0)
@m iy(n,a) ijy(n,a,1)
@m psi(n,a) psid(n,a,0)
@m dpsi(n,a) psid(n,a,1)
@ Subroutine to calculate the potentials. This is a straight
implementation of equations (31) and (22).
@<Functions...@>=
subroutine potential(ua,du,c,fm,chil,nmax,nmax1,legjy,psid,fla,ijy)
implicit_none
integer p0,p02,p022,p1,p11
parameter (p0=0, p02=1, p022=2, p1=3, p11=4)
integer j0,j1,j2,j02,j11,j22,j022
parameter (j0=0, j1=1, j2=2, j02=3, j11=4, j22=5, j022=6)
integer nmax,nmax1,n,a
real ua,du,c,fm,chil,legjy
dimension ua(0:nmax-1),fm(0:nmax-1),chil(0:nmax-1),
$ legjy(0:nmax1-1,0:6,0:3)
real psid
dimension psid(0:nmax1-1,0:4,0:1)
real fla,ijy
dimension fla(0:nmax-1),ijy(0:nmax1,0:6,0:1)
@<Define various indefinite integrals@>
do n=0,nmax-1
psi(n,p0) = y(n,j0)*ij(n,j0) + iy(n,j0)*j(n,j0)
psi(n,p02) = y(n,j0)*ij(n,j02) + y(n,j02)*ij(n,j2) +
$ iy(n,j0)*j(n,j02) + iy(n,j02)*j(n,j2)
psi(n,p022) = y(n,j0)*ij(n,j022) + y(n,j02)*ij(n,j22) +
$ y(n,j022)*ij(n,j2) + @|
$ iy(n,j0)*j(n,j022) + iy(n,j02)*j(n,j22) + iy(n,j022)*j(n,j2)
psi(n,p1) = y(n,j1)*ij(n,j1) + iy(n,j1)*j(n,j1)
psi(n,p11) = y(n,j1)*ij(n,j11) + y(n,j11)*ij(n,j1) +
$ iy(n,j1)*j(n,j11) + iy(n,j11)*j(n,j1)
dpsi(n,p0) = dy(n,j0)*ij(n,j0) + iy(n,j0)*dj(n,j0)
dpsi(n,p02) = dy(n,j0)*ij(n,j02) + dy(n,j02)*ij(n,j2) +
$ iy(n,j0)*dj(n,j02) + iy(n,j02)*dj(n,j2)
dpsi(n,p022) = dy(n,j0)*ij(n,j022) + dy(n,j02)*ij(n,j22) +
$ dy(n,j022)*ij(n,j2) + @|
$ iy(n,j0)*dj(n,j022) + iy(n,j02)*dj(n,j22) + iy(n,j022)*dj(n,j2)
dpsi(n,p1) = dy(n,j1)*ij(n,j1) + iy(n,j1)*dj(n,j1)
dpsi(n,p11) = dy(n,j1)*ij(n,j11) + dy(n,j11)*ij(n,j1) +
$ iy(n,j1)*dj(n,j11) + iy(n,j11)*dj(n,j1)
end do
return
end
@ Workspace for |potential| (|fla| and |ijy|).
@<Index spec...@>=
workalloc(nmax + (nmax1+1)*7*2)
@ Define various indefinite integrals needed in the definition of the
potentials. First do the integration to cell edges and then interpolate to
the cell centers.
@<Define various indefinite integrals@>=
do n=0,nmax-1
fla(n)=fm(n)*chil(n)*ua(n)^2/sqrt(1+(ua(n)/c)^2)*du/2
end do
do a=0,6
do n=0,nmax-1 /* Define the integrand */
ij(n+1,a)=fla(n)*j(n,a)
iy(n,a)=fla(n)*y(n,a)
end do
ij(0,a)=0 /* Define the boundary condtions */
iy(nmax,a)=0
end do
do n=1,nmax /* Integrate to cell edges */
do a=0,6
ij(n,a)=ij(n-1,a)+ij(n,a)
end do
end do
do n=nmax-1,0,-1
do a=0,6
iy(n,a)=iy(n+1,a)+iy(n,a)
end do
end do
do a=0,6
do n=0,nmax-1
ij(n,a)=ij(n,a)+ij(n+1,a)
iy(n,a)=iy(n,a)+iy(n+1,a)
end do
end do
@* Legendre harmonic decomposition. Define arrays for Legendre harmonics
and decompostion of $\chi$.
@<Glob...@>=
real chil,legp /* Legendre decomposition of $\chi$ */
dimension chil(0:nmax1-1,1:lmaxa1),legp(0:imax1-1,1:lmaxa1)
external initp
@
@<Index decl...@>=
integer ptr(chil),ptr(legp)
@
@<Index spec...@>=
alloc(chil,nmax1*lmaxa1)
alloc(legp,imax1*lmaxa1)
@ Initialize |legp|.
@<Init...@>=
call initp(csth0a,imax,imax1,lmaxa,lmaxa1,legp,work)
@ Define |initp|. Use the recursion formula for Legendre harmonics
$$ (l+1) P_{l+1}(\mu) = (2l+1)\mu P_l(\mu) - l P_{l-1}(\mu). $$
@<Functions...@>=
subroutine initp(x,imax,imax1,lmaxa,lmaxa1,legp,leg0)
implicit_none
integer imax,imax1,lmaxa,lmaxa1,i,l,la
real x,legp
dimension x(0:imax-1),legp(0:imax1-1,1:lmaxa1)
real leg0
dimension leg0(0:imax-1)
if (lmaxa<1) return /* Nothing to do */
la=1
do i=0,imax-1
leg0(i)=1
legp(i,la)=x(i)
end do
do la=2,lmaxa
l=lf(la)
do i=0,imax-1
leg0(i)=((2*l-3)*x(i)*legp(i,la-1)-(l-2)*leg0(i))/(l-1)
/* $P_{l-1}$ */
legp(i,la)=((2*l-1)*x(i)*leg0(i)-(l-1)*legp(i,la-1))/l
/* $P_{l}$ */
end do
end do
return
end
@ Workspace for |initp| (|leg0|).
@<Index spec...@>=
workalloc(imax)
@ Decompose $\chi$ into Legendre components. We have
$$
\chi_l(u)=\frac{2l+1}{2}\int_0^{\pi} \chi(u,\theta_0) P_l(\cos\theta_0)
\sin\theta_0\, d\theta_0.
$$
Since $\chi$ is odd in $\theta_0$ and it vanishes for the trapped
population this becomes
$$
\chi_l(u)=(2l+1)\int_0^{\theta_{0,\rm max}} \chi(u,\theta_0) P_l(\cos\theta_0)
\sin\theta_0\, d\theta_0.
$$
@<Functions...@>=
subroutine decomp(chi,nmax,nmax1,imax,imax1,lmaxa,lmaxa1,
$ th0a,dth0,legp,chil)
implicit_none
integer nmax,nmax1,imax,imax1,lmaxa,lmaxa1,la,l,n,i
real chi,th0a,dth0,legp,chil,sn
dimension chi(0:imax1-1,0:nmax-1),th0a(0:imax-1),
$ legp(0:imax1-1,1:lmaxa1),chil(0:nmax1-1,1:lmaxa1)
do la=1,lmaxa
l=lf(la)
do n=0,nmax-1
chil(n,la)=0
end do
do i=0,imax-1
sn=sin(th0a(i))*legp(i,la)
do n=0,nmax-1
chil(n,la)=chil(n,la)+chi(i,n)*sn
end do
end do
do n=0,nmax-1
chil(n,la)=(2*l+1)*dth0*chil(n,la)
end do
end do
return
end
@* Coefficients for isotropic collision operator. I.e., collisions off a
Maxwellian background.
@<Glob...@>=
real t,ze,zi /* Plasma parameters */
real duug,fu,dtt
dimension duug(0:nmax),fu(0:nmax-1),dtt(0:nmax-1)
real fm /* Maxwellian */
dimension fm(0:nmax-1)
real psid /* Arrays for potentials and their derivatives */
dimension psid(0:nmax1-1,0:4,0:1)
external maxwellian,isotropic
@
@<Index decl...@>=
integer ptr(duug),ptr(fu),ptr(dtt),ptr(fm),ptr(psid)
@
@<Index spec...@>=
alloc(duug,nmax+1)
alloc(fu,nmax)
alloc(dtt,nmax)
alloc(fm,nmax)
alloc(psid,nmax1*5*2)
@ Initialize isotropic terms.
@<Init...@>=
call maxwellian(ua,du,c,nmax,fm,t)
do n=0,nmax-1
chil(n,1)=1
end do
call isotropic(ua,du,c,fm,chil(0,1),nmax,nmax1,legjy(0,0,0,0),psid,
$ ze,zi,duug,fu,dtt,work)
@ This function calculates $D_{uu}$, $D_{\theta\theta}$ and $F_u$ from the
potentials. This is a transcription of equations (33).
@<Functions...@>=
subroutine isotropic(ua,du,c,fm,chil,nmax,nmax1,legjy,psid,
$ ze,zi,duug,fu,dtt,work)
implicit_none
integer p0,p02,p022,p1,p11
parameter (p0=0, p02=1, p022=2, p1=3, p11=4)
integer nmax,nmax1,n
real ua,c,psid,duug,fu,dtt,pi,g,u,ze,zi,du,fm,chil,legjy
dimension ua(0:nmax-1),fm(0:nmax-1),chil(0:nmax-1),
$ legjy(0:nmax1-1,0:6,0:3),psid(0:nmax1-1,0:4,0:1),
$ duug(0:nmax),fu(0:nmax-1),dtt(0:nmax-1)
real work
dimension work(0:nmax+(nmax1+1)*7*2-1)
external potential
call potential(ua,du,c,fm,chil,nmax,nmax1,legjy,psid,work(0),work(nmax))
pi=4*atan(const(1.0))
do n=0,nmax-1
u=ua(n)
g=sqrt(1+(u/c)^2)
duug(n)=ze * 4*pi*g/u*(2*g^2*dpsi(n,p02)-u*psi(n,p0)- @|
$ 8*g^2/c^2*dpsi(n,p022)+8*u/c^4*psi(n,p022))
dtt(n)=ze * 4*pi/(g*u)*(-g^2*dpsi(n,p02)-u/c^2*psi(n,p02)+ @|
$ 4*g^2/c^2*dpsi(n,p022)-4*u/c^4*psi(n,p022))
$ + zi * g/(2*u)
fu(n)=ze * 4*pi*g*(-dpsi(n,p1)+2/c^2*dpsi(n,p11))
end do
/* Interpolate |duug| to cell edges. Assume a constant first derivative at
|umax| and a constant value at $u=0$. */
duug(nmax)=(3*duug(nmax-1)-duug(nmax-2))/2
do n=nmax-1,1,-1
duug(n)=(duug(n)+duug(n-1))/2
end do
duug(0)=duug(0)
return
end
@ Define a Maxwellian distribution
$$
f_m(u) \propto \exp\biggl(-\frac{mc^2\gamma}{T}\biggr).
$$
@<Functions...@>=
subroutine maxwellian(ua,du,c,nmax,fm,t)
implicit_none
integer nmax,n
real ua,du,c,fm,t,sum,pi
dimension ua(0:nmax-1),fm(0:nmax-1)
sum=0
do n=0,nmax-1
fm(n)=exp(-ua(n)^2/(sqrt(1+(ua(n)/c)^2)+1)/t)
sum=sum+ua(n)^2*fm(n)
end do
pi=4*atan(const(1.0))
sum=4*pi*du*sum
do n=0,nmax-1
fm(n)=fm(n)/sum /* Normalize to unit density */
end do
return
end
@* Reaction terms. Calculate the term
$$
I(\chi_l(u))=
\frac{C^{e/e}(f_{m}(u),f_{m}(u)\chi_{l}(u)P_l(\cos\theta))}
{f_{m}(u)P_l(\cos\theta)},
$$
using equation (37). In fact we use a generalization of this equation for
arbitrary $l$. This term is added to the |reactl|.
@<Functions...@>=
subroutine reaction(ua,du,c,fm,t,chil,nmax,nmax1,legjy,psid,ze,reactl,l,
$ work)
implicit_none
integer p0,p02,p022,p1,p11
parameter (p0=0, p02=1, p022=2, p1=3, p11=4)
integer nmax,nmax1,n,l
real ua,du,c,legjy,psid,ze,reactl,fm,chil,t,pi,u,g
dimension ua(0:nmax-1),legjy(0:nmax1-1,0:6,0:3),psid(0:nmax1-1,0:4,0:1),
$ reactl(0:nmax-1),fm(0:nmax-1),chil(0:nmax-1)
real work
dimension work(0:nmax+(nmax1+1)*7*2-1)
external potential
call potential(ua,du,c,fm,chil,nmax,nmax1,legjy,psid,work(0),work(nmax))
pi=4*atan(const(1.0))
do n=0,nmax-1
u=ua(n)
g=sqrt(1+(u/c)^2)
reactl(n)=reactl(n)+ze*4*pi*(
$ fm(n)*chil(n)/g- @|
$ u/t*dpsi(n,p1)-2/(c^2*g)*psi(n,p1)+2*u/(c^2*t)*dpsi(n,p11)+ @|
$ u/t*dpsi(n,p0)-(u^2/(g*t^2)-1/t)*psi(n,p0)+ @|
$ (2*g*u/t^2-2*u/(c^2*t))*dpsi(n,p02)-
$ (l*(l+1)/(g*t^2)-2/(c^2*t))*psi(n,p02)- @|
$ 8*g*u/(c^2*t^2)*dpsi(n,p022)+
$ 4*((l+2)*(l-1)/g+2*g)/(c^2*t^2)*psi(n,p022)
$ )
end do
return
end
@ Sum up reaction terms
$$ \frac1{\lambda_1}
\sum_l^{l_{\rm max}}\sum_{l',l''\ \rm odd}^l
H_{l,l',l''} P_{l',0} I_l(\chi_{l'',0}).
$$
We can rewrite this as
$$ \frac1{\lambda_1}
\sum_{l'=1}^{l_{\rm max}} P_{l',0}
\sum_{l=l'}^{l_{\rm max}}
I_l\biggl( \sum_{l''=1}^l H_{l,l',l''} \chi_{l'',0} \biggr).
$$
@<Functions...@>=
subroutine reactall(ua,du,c,fm,t,chil,nmax,nmax1,imax,imax1,lam1a,
$ legjy,psid,legp,harr,ze,lmaxa,lmaxa1,react,reactl,chila,work)
implicit_none
integer nmax,nmax1,imax,imax1,lmaxa,lmaxa1
real ua,du,c,fm,t,chil,lam1a,legjy,psid,legp,harr,ze,react
dimension ua(0:nmax-1),fm(0:nmax-1),chil(0:nmax1-1,1:lmaxa1),
$ lam1a(0:imax-1),legjy(0:nmax1-1,0:6,0:3,1:lmaxa1),
$ psid(0:nmax1-1,0:4,0:1),legp(0:imax1-1,1:lmaxa1),
$ hv(lmaxa1,lmaxa1,lmaxa1),react(0:imax1-1,0:nmax-1)
integer la,la1,la2,n,i
real reactl,chila,work
dimension reactl(0:nmax-1),chila(0:nmax-1),work(0:nmax+(nmax1+1)*7*2-1)
external reaction
do n=0,nmax-1
do i=0,imax-1
react(i,n)=0
end do
end do
do la1=1,lmaxa
@<Compute innermost sums@>
do n=0,nmax-1
do i=0,imax-1
react(i,n)=react(i,n)+reactl(n)*legp(i,la1)
end do
end do
end do
do n=0,nmax-1
do i=0,imax-1
react(i,n)=react(i,n)/lam1a(i)
end do
end do
return
end
@ Workspace for |reactall| (|reactl| and |chila|) and for |potential|
(|fla| and |ijy|).
@<Index spec...@>=
workalloc(nmax * 2 + nmax + (nmax1+1)*7*2)
@ Do the innermost two sums.
@<Compute innermost sums@>=
do n=0,nmax-1
reactl(n)=0
end do
do la=la1,lmaxa
do n=0,nmax-1
chila(n)=0
end do
do la2=1,la
do n=0,nmax-1
chila(n)=chila(n)+hv(la,la1,la2)*chil(n,la2)
end do
end do
call reaction(ua,du,c,fm,t,chila,nmax,nmax1,
$ legjy(0,0,0,la),psid,ze,reactl,lf(la),work)
end do
@* Inhomogeneous magnetic field.
@<Glob...@>=
real eps,deltab /* Inverse aspect ratio */
external deltab /* Variation of magnetic field */
@ Define magnetic field variation,
$$b=\frac{B(\phi)}{B(0)},\qquad \delta b=b-1,$$
where $\phi$ is the poloidal angle. We assume that
$$b=\frac{1+\epsilon}{1+\epsilon\cos\phi};\qquad
\delta b=\frac{2\epsilon\sin^2(\phi/2)}{1+\epsilon\cos\phi}.$$
@<Functions...@>=
function deltab(phi,eps)
implicit_none
real deltab,phi,eps
deltab=2*eps*sin(phi/2)^2/(1+eps*cos(phi))
return
end
@ Arrays to hold various averaged quantities.
@<Glob...@>=
integer kmax
real lam1,lam2,lam1a,lam2g,harr
dimension lam1a(0:imax-1),lam2g(0:imax),hv(lmaxa1,lmaxa1,lmaxa1)
external lam1,lam2,hinit
@
@<Index decl...@>=
integer ptr(lam1a),ptr(lam2g),ptr(harr)
@
@<Index spec...@>=
alloc(lam1a,imax)
alloc(lam2g,imax+1)
alloc(harr,hind(lmaxa1,lmaxa1,lmaxa1))
@ Initialize these arrays.
@<Init...@>=
do i=0,imax-1
lam1a(i)=lam1(eps,th0a(i),kmax)
end do
do i=0,imax
lam2g(i)=lam2(eps,th0g(i),kmax)
end do
call hinit(eps,kmax,lmaxa1,harr,work)
@ Define averages over the magnetic field. We need three:
$$\eqalign{
\lambda_1&=\frac 1{2\pi} \int d\phi \frac{\cos\theta_0}{\cos\theta},\cr
\lambda_2&=\frac 1{2\pi}
\int \frac{d\phi}{b} \frac{\cos\theta}{\cos\theta_0},\cr
H_{l,l',l''}&=\frac 1{2\pi} \int d\phi\,b G_{l,l'} G_{l,l''},\cr
}$$
where
$$\sin^2\theta=b \sin^2\theta_0,$$
and hence
$$\cos\theta=\cos\theta_0\sqrt{1-\delta b\tan^2\theta_0}.$$
The following functions are specialized to the case of passing particles
only.
@ Define $\lambda_1$.
@<Functions...@>=
function lam1(eps,th0,kmax)
implicit_none
integer kmax,k
real lam1,eps,th0,phi,deltab,deltabv,sum,costh,pi
external deltab
pi=4*atan(const(1.0))
sum=0
do k=0,kmax-1
phi=2*pi*(k+const(0.25))/kmax
deltabv=deltab(phi,eps)
costh=sqrt(1-min(const(1.0),deltabv*tan(th0)^2))
sum=sum+1/costh
end do
lam1=sum/kmax
return
end
@ Define $\lambda_2$.
@<Functions...@>=
function lam2(eps,th0,kmax)
implicit_none
integer kmax,k
real lam2,eps,th0,phi,deltab,deltabv,sum,costh,pi
external deltab
pi=4*atan(const(1.0))
sum=0
do k=0,kmax-1
phi=2*pi*(k+const(0.25))/kmax
deltabv=deltab(phi,eps)
costh=sqrt(1-min(const(1.0),deltabv*tan(th0)^2))
sum=sum+costh/(1+deltabv)
end do
lam2=sum/kmax
return
end
@ Define $H_{l,l',l''}$ via
$$
H_{l,l',l''}=\frac{2l+1}{2l''+1}
\lave{ b G_{l,l'} G_{l,l''}}.
$$
|gv| allows us to index into |garr| conveniently. Enough space is allowed
to accomodate |gv(la,0)|, |gv(la,la+1)|, and |gv(la,la+2)|. These can be
initialized to $0$ thereby simplifying the recursion for $G$.
$H$ is defined out to |lmaxa1| because we need $H_{1,1,1}$ even when
$l_{\rm max}=0$.
@<Functions...@>=
subroutine hinit(eps,kmax,lmaxa1,harr,garr)
implicit_none
integer kmax,k,lmaxa1,la,la1,la2,l,l1
real harr,eps,phi,deltab,b,pi,garr
dimension hv(lmaxa1,lmaxa1,lmaxa1)
dimension gv(lmaxa1,lmaxa1)
external deltab
pi=4*atan(const(1.0))
@<Set initialial values for |harr| and |garr|@>
do k=0,kmax-1
phi=2*pi*(k+const(0.25))/kmax
b=1+deltab(phi,eps)
@<Fill |garr|@>
do la=1,lmaxa1
do la1=1,la
do la2=1,la
hv(la,la1,la2)=hv(la,la1,la2)+b*gv(la,la1)*gv(la,la2)
end do
end do
end do
end do
do la=1,lmaxa1
do la1=1,la
do la2=1,la
hv(la,la1,la2)=(2*lf(la)+1)/areal(2*lf(la2)+1)*
$ hv(la,la1,la2)/kmax
end do
end do
end do
return
end
@ Workspace for |hinit| (|garr|).
@<Index spec...@>=
workalloc(gind(lmaxa1,lmaxa1))
@ Initialize |harr| and set the edge values of |garr| to zero.
@<Set initialial values for |harr| and |garr|@>=
do la=1,lmaxa1
do la1=1,la
do la2=1,la
hv(la,la1,la2)=0
end do
end do
end do
do la=-2,lmaxa1-1
gv(la,0)=0
end do
do la=-1,lmaxa1-1
gv(la,la+1)=0
gv(la,la+2)=0
end do
gv(1,1)=1
@ Use the recursion relation for $G$ to fill |garr|.
@<Fill |garr|@>=
do la=2,lmaxa1
l=lf(la)
do la1=1,la
l1=lf(la1)
gv(la,la1)=(2*(2*l-3)*(l^2-3*l+1)*gv(la-1,la1) - @|
$ (l-3)*(l-2)*(2*l-1)*gv(la-2,la1)) / (l*(l-1)*(2*l-5)) + @|
$ (2*l-3)*(2*l-1)*b*
$ ( (l1+1)*(l1+2) * gv(la-1,la1+1)/ ((2*l1+3)*(2*l1+5)) - @|
$ 2*(l1^2+l1-1) * gv(la-1,la1)/((2*l1-1)*(2*l1+3)) + @|
$ l1*(l1-1) * gv(la-1,la1-1)/((2*l1-3)*(2*l1-1)) )/(l*(l-1))
end do
end do
@* Residue calculation. The array |resid| is already initialized with the
reaction term. Two methods for the calculation are given here: the
straightforward coding of the residue term, which is slow but easy to
check; and a faster method using the arrays |uinv| and |th0inv| which are
used in the tridiagonal elimination.
@m RESIDUE 0 /* Should we calculate the residue in the straightforward way? */
@<Functions...@>=
@#if RESIDUE
subroutine residue(du,ua,ug,nmax,dth0,th0a,th0g,csth0a,imax,imax1,
$ chi,duug,fu,dtt,c,lam1a,lam2g,rho,resid)
implicit_none
integer nmax,imax,imax1,n,i
real du,ua,ug,dth0,th0a,th0g,csth0a,chi,duug,fu,dtt,c,
$ lam1a,lam2g,rho,resid
dimension ua(0:nmax-1),ug(0:nmax),th0a(0:imax-1),th0g(0:imax),
$ csth0a(0:imax-1),
$ chi(0:imax1-1,0:nmax-1),duug(0:nmax),fu(0:nmax-1),dtt(0:nmax-1),
$ lam1a(0:imax-1),lam2g(0:imax),resid(0:imax1-1,0:nmax-1)
@<Energy diffusion term@>
@<Friction term@>
@<Pitch-angle scattering term@>
@<Electric field term@>
/* Normalize to pitch-angle diffusion coefficient. */
do n=0,nmax-1
do i=0,imax-1
resid(i,n)=resid(i,n)*ua(n)^2/dtt(n)
end do
end do
return
end
@#endif
@ The faster way.
@<Functions...@>=
@#if !RESIDUE
@#if VAX
options /check=nobounds
@#endif
subroutine residue(ua,nmax,nmax1,csth0a,imax,imax1,chi,dtt,c,lam1a,
$ uinv,th0inv,resid)
implicit_none
integer nmax,nmax1,imax,imax1,n,i
real ua,csth0a,chi,dtt,c,lam1a,resid,uinv,th0inv
dimension ua(0:nmax-1),chi(0:imax1-1,0:nmax-1),csth0a(0:imax-1),
$ dtt(0:nmax-1),lam1a(0:imax-1),resid(0:imax1-1,0:nmax-1),
$ uinv(0:nmax1-1,-1:1),th0inv(0:imax1-1,-1:1)
@<Electric field term@>
/* Normalize to pitch-angle diffusion coefficient and add in
differential terms. */
do n=0,nmax-1
do i=0,imax-1
resid(i,n)=resid(i,n)*ua(n)^2/dtt(n)+ @|
$ chi(i,n-1)*uinv(n,-1)+chi(i,n)*uinv(n,0)+
$ chi(i,n+1)*uinv(n,1)+ @|
$ chi(i-1,n)*th0inv(i,-1)+chi(i,n)*th0inv(i,0)+
$ chi(i+1,n)*th0inv(i,1)
end do
end do
return
end
@#endif
@ Calculate contributions to the residue from the energy diffusion term
$$\frac{1}{u^2}\frac{\d}{\d u} u^2 D_{uu} \frac{\d \chi}{\d u}.$$
Boundary at $u=0$ takes care of itself. Boundary at $u=u_{\rm max}$
handled by specifying
|chi(i,nmax)=chi(i,nmax-1)+rho*(chi(i,nmax-1)-chi(i,nmax-2))|. For
$\rho=1$, we have $\chi''(u_{\rm max})=0$, and for $\rho=0$, we have
$\chi'(u_{\rm max})=0$.
@<Energy diffusion term@>=
do i=0,imax-1
n=0
resid(i,n)=resid(i,n) + @|
$ (ug(n+1)^2*duug(n+1)*(chi(i,n+1)-chi(i,n))/du)/(ua(n)^2*du)
do n=1,nmax-2
resid(i,n)=resid(i,n) + @|
$ (ug(n+1)^2*duug(n+1)*(chi(i,n+1)-chi(i,n))/du - @|
$ ug(n)^2*duug(n)*(chi(i,n)-chi(i,n-1))/du)/(ua(n)^2*du)
end do
n=nmax-1
resid(i,n)=resid(i,n) + @|
$ (ug(n+1)^2*duug(n+1)*rho*(chi(i,n)-chi(i,n-1))/du - @|
$ ug(n)^2*duug(n)*(chi(i,n)-chi(i,n-1))/du)/(ua(n)^2*du)
end do
@ Calculate contributions to the residue from the friction term
$$F_{u} \frac{\d \chi}{\d u}.$$
Boundary at $u=0$ handled by assuming $\chi$ is odd in $u$, i.e.,
|chi(i,-1)=-chi(i,0)|. Boundary at $u=u_{\rm max}$ handled in the same way
as in the calculation of the energy diffusion term.
@<Friction term@>=
do i=0,imax-1
n=0
resid(i,n)=resid(i,n)+fu(n)*(chi(i,n+1)+chi(i,n))/(2*du)
do n=1,nmax-2
resid(i,n)=resid(i,n)+fu(n)*(chi(i,n+1)-chi(i,n-1))/(2*du)
end do
n=nmax-1
resid(i,n)=resid(i,n)+fu(n)*(1+rho)*(chi(i,n)-chi(i,n-1))/(2*du)
end do
@ Calculate contributions to the residue from the pitch-angle scattering
term
$$\frac{D_{\theta\theta}}{u^2}
\frac{1}{\lambda_1\sin\theta_0}\frac{\d}{\d\theta_0} \sin\theta_0 \lambda_2
\frac{\d}{\d\theta_0}.$$
Boundary at $\theta_0=0$ takes care of itself. At $\theta_{0,\rm max}
=\theta_t$, assume that $\chi=0$.
@<Pitch-angle scattering term@>=
do n=0,nmax-1
i=0
resid(i,n)=resid(i,n)+dtt(n)/ua(n)^2 * @|
$ (sin(th0g(i+1))*lam2g(i+1)*(chi(i+1,n)-chi(i,n))/dth0)/
$ (lam1a(i)*sin(th0a(i))*dth0)
do i=1,imax-2
resid(i,n)=resid(i,n)+dtt(n)/ua(n)^2 * @|
$ (sin(th0g(i+1))*lam2g(i+1)*(chi(i+1,n)-chi(i,n))/dth0 - @|
$ sin(th0g(i))*lam2g(i)*(chi(i,n)-chi(i-1,n))/dth0)/
$ (lam1a(i)*sin(th0a(i))*dth0)
end do
i=imax-1
resid(i,n)=resid(i,n)+dtt(n)/ua(n)^2 * @|
$ (sin(th0g(i+1))*lam2g(i+1)*(-chi(i,n))/(dth0/2) - @|
$ sin(th0g(i))*lam2g(i)*(chi(i,n)-chi(i-1,n))/dth0)/
$ (lam1a(i)*sin(th0a(i))*dth0)
end do
@ Finally add in the electric field term:
$$\frac{v \cos\theta_0}{\lambda_1}. $$
@<Electric field term@>=
do n=0,nmax-1
do i=0,imax-1
resid(i,n)=resid(i,n)+
$ csth0a(i)/lam1a(i)*ua(n)/sqrt(1+(ua(n)/c)^2)
end do
end do
@ Calculate norm of residue.
@<Functions...@>=
subroutine error(du,ua,umax,nmax,dth0,th0a,th0max,imax,imax1,resid,chi,
$ abserr,relerr,aerr,rerr)
implicit_none
integer nmax,imax,imax1,i,n
real du,ua,umax,dth0,th0a,th0max,resid,chi,abserr,relerr
dimension ua(0:nmax-1),th0a(0:imax-1),resid(0:imax1-1,0:nmax-1),
$ chi(0:imax1-1,0:nmax-1)
real aerr,rerr
dimension aerr(0:imax-1),rerr(0:imax-1)
do i=0,imax-1
aerr(i)=0
rerr(i)=0
end do
do n=0,nmax-1
do i=0,imax-1
aerr(i)=aerr(i)+resid(i,n)^2*ua(n)^2
rerr(i)=rerr(i)+(resid(i,n)/max(epsilon,abs(chi(i,n))))^2*ua(n)^2
end do
end do
abserr=0
relerr=0
do i=0,imax-1
abserr=abserr+sin(th0a(i))*aerr(i)
relerr=relerr+sin(th0a(i))*rerr(i)
end do
/* Normalize to total volume of velocity space */
abserr=sqrt(abserr*du*dth0/(umax^3/3*(1-cos(th0max))))
relerr=sqrt(relerr*du*dth0/(umax^3/3*(1-cos(th0max))))
return
end
@ Workspace for |error| (|aerr| and |rerr|).
@<Index spec...@>=
workalloc(2*imax)
@* Matrices for inversion.
@<Glob...@>=
real uinv,th0inv,rho
dimension uinv(0:nmax1-1,-1:1),th0inv(0:imax1-1,-1:1)
external matrixinit
@
@<Index decl...@>=
integer ptr(uinv),ptr(th0inv)
@
@<Index spec...@>=
alloc(uinv,nmax1*3)
alloc(th0inv,imax1*3)
@ Initialize matrices.
@<Init...@>=
call matrixinit(du,ua,ug,nmax,nmax1,dth0,th0a,th0g,imax,imax1,
$ duug,fu,dtt,lam1a,lam2g,uinv,th0inv,rho)
@ Calculate inversion matrices. These have to be compatible with
subroutine |residue|.
@<Functions...@>=
subroutine matrixinit(du,ua,ug,nmax,nmax1,dth0,th0a,th0g,imax,imax1,
$ duug,fu,dtt,lam1a,lam2g,uinv,th0inv,rho)
implicit_none
integer nmax,nmax1,imax,imax1,n,i
real du,ua,ug,dth0,th0a,th0g,duug,fu,dtt,lam1a,lam2g,
$ uinv,th0inv,rho
dimension ua(0:nmax-1),ug(0:nmax),th0a(0:imax-1),th0g(0:imax),
$ duug(0:nmax),fu(0:nmax-1),dtt(0:nmax-1),
$ lam1a(0:imax-1),lam2g(0:imax),
$ uinv(0:nmax1-1,-1:1),th0inv(0:imax1-1,-1:1)
@<Matrices for $u$ inversion@>
@<Matrices for $\theta_0$ inversion@>
return
end
@ Define matrices for inversion of diffusion and friction operator. These
are one-dimensional because the operator is independent of $\theta_0$.
@<Matrices for $u$ inversion@>=
do n=0,nmax-1
uinv(n,-1)=ug(n)^2*duug(n)/(ua(n)^2*du^2)
uinv(n,1)=ug(n+1)^2*duug(n+1)/(ua(n)^2*du^2)
uinv(n,0)=-uinv(n,-1)-uinv(n,1)
uinv(n,-1)=uinv(n,-1)-fu(n)/(2*du)
uinv(n,1)=uinv(n,1)+fu(n)/(2*du)
end do
n=0 /* Odd at $u=0$ */
uinv(n,0)=uinv(n,0)-uinv(n,-1)
uinv(n,-1)=0
n=nmax-1 /* Apply boundary condition at $u_{\rm max}$ */
uinv(n,-1)=uinv(n,-1)-rho*uinv(n,1)
uinv(n,0)=uinv(n,0)+(1+rho)*uinv(n,1)
uinv(n,1)=0
do n=0,nmax-1
uinv(n,-1)=uinv(n,-1)*ua(n)^2/dtt(n)
uinv(n,0)=uinv(n,0)*ua(n)^2/dtt(n)
uinv(n,1)=uinv(n,1)*ua(n)^2/dtt(n)
end do
@ Define matrices for inversion of pitch-angle scattering operator. These
are one-dimensional because, after divided out $D_{\theta\theta,0}/u^2$,
the operator is independent of $u$.
@<Matrices for $\theta_0$ inversion@>=
do i=0,imax-1
th0inv(i,-1)=sin(th0g(i))*lam2g(i)/(lam1a(i)*sin(th0a(i))*dth0^2)
th0inv(i,1)=sin(th0g(i+1))*lam2g(i+1)/(lam1a(i)*sin(th0a(i))*dth0^2)
th0inv(i,0)=-th0inv(i,-1)-th0inv(i,1)
end do
i=0 /* Nothing special needs to be done */
i=imax-1 /* Zero on trapped passing boundary */
th0inv(i,0)=th0inv(i,0)-th0inv(i,1)
th0inv(i,1)=0
@ Solve the tridiagonal equations.
This subroutine solves the tridiagonal equation
$$
x_i-\alpha\Delta t(a_i x_{i-1}+b_i x_i+c_i x_{i+1})=z_i,
$$
for $0\le i< n$, with $a_0=c_{n-1}=0$.
Vectorized to handle $m$ parallel systems. |ms| is the skipping distance
in the vectorizing direction. |ns| is the skipping distance in the
solution direction between $x_i$ and $x_{i+1}$. In this version, |a|, |b|, and
|c| are the same for each system.
Avoid using any working storage by clobbering |z| and using |x| as a
temporary.
@<Functions...@>=
@#if VAX
options /check=nobounds
@#endif
subroutine solve(x,ns,n,ms,m,a,b,c,z,dt,alpha)
implicit_none
integer n,ns,m,ms,i,k
real x,a,b,c,z,dt,alpha,dt2,den
dimension x(0:ms-1,0:*),z(0:ms-1,0:*),a(0:n-1),b(0:n-1),c(0:n-1)
dt2=-alpha*dt
i=0
do k=0,m-1
den=1+dt2*b(i)
x(ns*i,k)=z(ns*i,k)/den
z(ns*i,k)=-dt2*c(i)/den
end do
do i=1,n-2
do k=0,m-1
den=1+dt2*(b(i)+a(i)*z(ns*(i-1),k))
x(ns*i,k)=(z(ns*i,k)-dt2*a(i)*x(ns*(i-1),k))/den
z(ns*i,k)=-dt2*c(i)/den
end do
end do
i=n-1
do k=0,m-1
den=1+dt2*(b(i)+a(i)*z(ns*(i-1),k))
x(ns*i,k)=(z(ns*i,k)-dt2*a(i)*x(ns*(i-1),k))/den
z(ns*i,k)=0
end do
/* The equations have now been converted to the form $x_i=r_i x_{i+1}+s_i$
(with $r=y$, $s=x$). */
do i=n-2,0,-1
do k=0,m-1
x(ns*i,k)=z(ns*i,k)*x(ns*(i+1),k)+x(ns*i,k)
end do
end do
return
end
@* Time advancement.
@<Glob...@>=
integer tmax,time,tskip
real resid,abserr,relerr,alpha,dt,resida
dimension resid(0:imax1-1,0:nmax-1),resida(0:imax1-1,0:nmax-1)
external decomp,reactall,residue,error,solve
@
@<Index decl...@>=
integer ptr(resid),ptr(resida)
@
@<Index spec...@>=
alloc(resid,imax1*nmax)
alloc(resida,imax1*nmax)
@ Calculate reaction term.
@<Make a step@>=
call decomp(chi,nmax,nmax1,imax,imax1,lmaxa,lmaxa1,th0a,dth0,legp,chil)
call reactall(ua,du,c,fm,t,chil,nmax,nmax1,imax,imax1,lam1a,
$ legjy(0,0,0,1),psid,legp,harr,ze,lmaxa,lmaxa1,resid,
$ work(0),work(nmax),work(2*nmax))
@ Calculate residue and its norm.
@<Make a step@>=
@#if RESIDUE
call residue(du,ua,ug,nmax,dth0,th0a,th0g,csth0a,imax,imax1,chi,
$ duug,fu,dtt,c,lam1a,lam2g,rho,resid)
@#else
call residue(ua,nmax,nmax1,csth0a,imax,imax1,chi,dtt,c,lam1a,
$ uinv,th0inv,resid)
@#endif
if (mod(time,tskip)==0) then
call error(du,ua,umax,nmax,dth0,th0a,th0max,imax,imax1,resid,chi,
$ abserr,relerr,work(0),work(imax))
end if
@ Advance $\chi$.
@<Make a step@>=
call solve(resida,imax1,nmax,1,imax,uinv(0,-1),uinv(0,0),uinv(0,1),
$ resid,dt,alpha)
call solve(resid,1,imax,imax1,nmax,th0inv(0,-1),th0inv(0,0),th0inv(0,1),
$ resida,dt,alpha)
do n=0,nmax-1
do i=0,imax-1
chi(i,n)=chi(i,n)+dt*resid(i,n)
end do
end do
@* Current and conductivity.
@<Glob...@>=
real curv,current,cond,e
external current
@ Compute current and conductivity. The conductivity is define as
$J_{0\parallel}/E_{0\parallel}$ where
$$E_{0\parallel} = B_0\frac{\int E_{\parallel} dl}{\int B\,dl}.$$
Substituting $E_\parallel = (B_\zeta/B) \Phi/(2\pi R)$ and $\Phi=TQ/L$ gives
$E_\parallel = T/\int b \,dl/L$.
@<Make a step@>=
if (mod(time,tskip)==0) then
curv=current(du,ua,c,nmax,dth0,th0a,imax,imax1,fm,chi,work)
e=t/hv(1,1,1)
cond=curv*zi/(t*sqrt(t)*e)
write(outunit,*) time,cond,abserr,relerr
end if
@ Function to calculate current.
@<Functions...@>=
function current(du,ua,c,nmax,dth0,th0a,imax,imax1,fm,chi,cur)
implicit_none
integer nmax,imax,imax1,n,i
real current,du,ua,c,dth0,th0a,fm,chi
dimension ua(0:nmax-1),th0a(0:imax-1),
$ fm(0:nmax-1),chi(0:imax1-1,0:nmax-1)
real cur,pi
dimension cur(0:imax-1)
pi=4*atan(const(1.0))
do i=0,imax-1
cur(i)=0
end do
do n=0,nmax-1
do i=0,imax-1
cur(i)=cur(i)+chi(i,n)*fm(n)*ua(n)**3/sqrt(1+(ua(n)/c)^2)
end do
end do
current=0
do i=0,imax-1
current=current+sin(th0a(i))*cos(th0a(i))*cur(i)
end do
current=4*pi*du*dth0*current
return
end
@ Workspace for |current| (|cur|).
@<Index spec...@>=
workalloc(imax)
@* Write out results to disk.
@<Write out results@>=
write(diskunit) nmax,imax,kmax,lmax
write(diskunit) single(t),single(c2),single(ze),single(zi),single(eps)
write(diskunit) single(umax),single(th0max),single(du),single(dth0),
$ single(rho)
write(diskunit) single(dt),tmax,single(alpha)
write(diskunit) single(cond),single(abserr),single(relerr)
write(diskunit) (single(ua(n)),n=0,nmax-1)
write(diskunit) (single(ug(n)),n=0,nmax)
write(diskunit) (single(th0a(n)),n=0,imax-1)
write(diskunit) (single(th0g(n)),n=0,imax)
write(diskunit) (single(fm(n)),n=0,nmax-1)
write(diskunit) ((single(chi(i,n)),i=0,imax-1),n=0,nmax-1)
@* Legendre functions. It is convenient to express the potentials in terms
of a set of functions defined by
$$
\eqalignno{
j\subb l1a(z) &=\sqrt{\frac{\pi}{2z}}P^{-l-1/2}_{a-1/2}(\sqrt{1+z^2}),\cr
y\subb l1a(z) &=(-1)^{l+1}\sqrt{\frac{\pi}{2z}}P^{l+1/2}_{a-1/2}(\sqrt{1+z^2}),
\cr}
$$
where $P^\mu_\nu$ is the associated Legendre function of the first kind.
Multiple index functions are defined recursively by
$$
j\subb l{k+2}{\s aa'} =
\cases{
\displaystyle\frac{j\subb l{k+1}{\s a}-j\subb l{k+1}{\s a'}}
{a^2-a'^2},&for $a\ne a'$,\cr
\displaystyle\pd{j\subb l{k+1}{\s a}}{(a^2)},
&for $a=a'$,\cr
}$$
with $y\subb lk\s$ defined in a similar fashion.
The two families of functions are related by
$$y\subb lk{\s}=(-1)^{l+1}j\subb {-l-1}k{\s}.$$
In fact we shall work with two derived functions, scaled Legendre functions
$\ja{l,a}(z)$ defined by
$$j\subb l1a(z) = \frac{z^l}{(2l+1)!!} \ja{l,a}(z),$$
and normalized Legendre functions defined by
$$ \jn lk{\s}(u) = c^{2k+l-2} j\subb lk{\s}(u/c). $$
For $k=1$, we have
$$\jn\subb l1a(z) = \frac{u^l}{(2l+1)!!} \ja{l,a}(u/c).$$
The methods used for computing these functions are given in the appendix.
@ Scaled Legendre functions. This subroutine returns the scaled Legendre
functions $\hbox{|jarray(l,a)|}=\ja{l,a}(z)$ for $-|lmax|-1 \le l \le
|lmax|$ and $0\le a \le |amax|$.
@<Functions...@>=
subroutine legendrejs(z,lmax,amax,jarray,lmax1)
implicit_none
integer l,a,lmax,amax,lmax1,lstart
real z,jarray,z2,g,ja0,ja1,ja2,ratio,sigma
dimension jarray(-lmax1-1:lmax1,0:amax)
if (amax<0 || lmax<0) return /* Nothing to do! */
if (lmax>lmax1) return /* Error: not enough room! */
z2=z^2
g=sqrt(1+z2)
@< Do the cases $-a\le l<a$ @>
@< Do the cases $l<-a$ @>
@< Do the cases $l\ge0$, $a=0$ @>
@< Do the other cases, i.e., $a>0$, $l\ge a$ @>
return
end
@ Calculate the scaled Legendre functions for the easy case, namely $-a\le
l<a$.
@< Do the cases $-a\le l<a$ @>=
do a=1,amax /* $a=0$ is excluded */
ja2=1 /* $\ja{a-1,a}$ */
ja1=g /* $\ja{a-2,a}$ */
do l=a-3,lmax-1,-1 @/
/* Get things started by iterating down to $|lmax|-1$ if need be.
This is from equation (A29) */
ja0=g*ja1-z2*(l-a+2)*(l+a+2)*ja2/((2*l+3)*(2*l+5))
ja2=ja1 /* $|ja2| = \ja{l+1,a}$ */
ja1=ja0 /* $|ja1| = \ja{l,a}$ */
end do
l=min(lmax,a-1)
jarray(l,a)=ja2
l=l-1
jarray(l,a)=ja1
do l=min(lmax,a-1)-2,max(-a,-lmax-1),-1 @/
/* Main iteration loop---again from equation (A29) */
jarray(l,a)=g*jarray(l+1,a)-
$ z2*(l-a+2)*(l+a+2)*jarray(l+2,a)/((2*l+3)*(2*l+5))
end do
end do
@ Calculate the scaled Legendre functions for the other easy case, namely
$l<-a$.
@< Do the cases $l<-a$ @>=
do a=0,min(amax,lmax)
l=-a-1
jarray(l,a)=1
l=l-1
if (l>=-lmax-1) then
jarray(l,a)=g
end if
do l=-a-3,-lmax-1,-1 @/
/* Main iteration loop---again from equation (A29) */
jarray(l,a)=g*jarray(l+1,a)-
$ z2*(l-a+2)*(l+a+2)*jarray(l+2,a)/((2*l+3)*(2*l+5))
end do
end do
@ Calculate $\ja{l\ge0,0}$. Depending on the magnitude of the argument we
either use backwards recursion (for $z \le |lmax|+1$) or forwards
recursion using $\ja{-1,0}=1$ and $\ja{0,0} = (\sinh^{-1} z)/z$.
@< Do the cases $l\ge0$, $a=0$ @>=
a = 0
if (abs(z) <= lmax+1) then @/
/* Estimate where to start the recursion. The 10 is a extra guard. */
lstart = lmax + log(epsilon)/(2*log((epsilon+abs(z))/(g+1))) + 10 @/
/* Initial guess for ratio $\ja{L'+1,a}/\ja{L',a}$. */
ratio = 2/(g+1)
do l=lstart,lmax-1,-1
ratio=1/(g - z2*(l-a+2)*(l+a+2)*ratio/((2*l+3)*(2*l+5))) @/
/* At this point, |ratio| is approximation for $\ja{l+1,a}/\ja{l,a}$. */
end do
jarray(lmax,a)=ratio
ja0=1
do l=lmax-2,-1,-1
jarray(l+1,a)=ja0
ja0=g*jarray(l+1,a)-
$ z2*(l-a+2)*(l+a+2)*jarray(l+2,a)/((2*l+3)*(2*l+5))
end do
ja0=1/ja0
do l=0,lmax
jarray(l,a)=ja0*jarray(l,a)
end do
else
/* Use forward recurrence for $z>|lmax|+1$. */
sigma=log(g+abs(z)) /* $\sinh^{-1}(\abs z)$ */
l=0
jarray(l,a)=sigma/abs(z) /* Starting value. We assume
that $\ja{-1,0}$ has already been set to $1$. */
do l=1,lmax @/
/* Equation (A29) */
jarray(l,a)=(g*jarray(l-1,a)-jarray(l-2,a))*
$ (2*l-1)*(2*l+1)/(z2*(l-a)*(l+a))
end do
end if
@ Do the other cases, i.e., $a>0$, $l\ge a$.
@< Do the other cases, i.e., $a>0$, $l\ge a$ @>=
do a=1,min(lmax,amax)
l=lmax @/
/* Initialize backwards recursion using equation (A30). */
jarray(l,a)=((2*l+1)*jarray(l-1,a-1)-
$ (l+1-a)*g*jarray(l,a-1))/(l+a)
l=lmax-1
if (l==0) then /* Implies that $a=1$. Use $\ja{-1,0} = 1$ */
ja0=((2*l+1)-(l+1-a)*g*jarray(l,a-1))/(l+a)
else
ja0=((2*l+1)*jarray(l-1,a-1)-
$ (l+1-a)*g*jarray(l,a-1))/(l+a)
end if
do l=lmax-2,a-1,-1
jarray(l+1,a)=ja0
ja0=g*jarray(l+1,a)-
$ z2*(l-a+2)*(l+a+2)*jarray(l+2,a)/((2*l+3)*(2*l+5))
end do
ja0=1/ja0
do l=a,lmax
jarray(l,a)=ja0*jarray(l,a)
end do
end do
@ Legendre functions of the first kind. This returns $
|jarray(l,a)|=j\subb l1a(z)$ for $0\le l \le |lmax|$ and $0\le a \le
|amax|$. In fact, we don't need this function; so it's commented out.
@<Functions...@>=
@#if 0
subroutine legendrej(z,lmax,amax,jarray,lmax1)
implicit_none
integer l,a,lmax,amax,lmax1
real z,jarray,scale
dimension jarray(-lmax1-1:lmax1,0:amax)
external legendrejs
call legendrejs(z,lmax,amax,jarray,lmax1)
scale=1
do l=1,lmax
scale=scale*z/(2*l+1)
do a=0,amax
jarray(l,a)=scale*jarray(l,a)
end do
end do
scale=1
do l=-1,-lmax-1,-1
scale=scale/z*(2*l+3)
do a=0,amax
jarray(l,a)=scale*jarray(l,a)
end do
end do
return
end
@#endif
@ Calculate normalized Legendre functions needed in computation of
collision operator. The single-index functions
$\jn l1a(u)$ are given in terms of $\ja{l,a}(u/c)$.
The multiple-index function are given in terms of the single-index
functions by by equations (A26):
$$\eqalign{
\jn l2{02}(u)&=\frac u2 \jn{l+1}11(u),\cr
\jn l2{11}(u)&=\frac u2 \jn{l+1}10(u),\cr
\jn l2{22}(u)&=\frac u2 \jn{l+1}11(u)+\frac{u^2}{2c^2} \jn{l+2}10(u),\cr
\jn l3{022}(u)&=\frac {u^2}8 \jn{l+2}10(u).\cr
}$$
The derivatives are given by equation (A21):
$$
\frac d{du}\jn lk*(u)=
\frac1{\gamma}\jn{l-1}k*(u)-\frac{l+1}{u}\jn lk*(u).
$$
On output, this returns
$$\eqalign{
\hbox{|jarray(l,0)|} &= \jn l10(u), \cr
\hbox{|jarray(l,1)|} &= \jn l11(u), \cr
\hbox{|jarray(l,2)|} &= \jn l12(u), \cr
\hbox{|jarray(l,3)|} &= \jn l2{02}(u), \cr
\hbox{|jarray(l,4)|} &= \jn l2{11}(u), \cr
\hbox{|jarray(l,5)|} &= \jn l2{22}(u), \cr
\hbox{|jarray(l,6)|} &= \jn l3{022}(u), \cr
\hbox{|djarray(l,a)|} &=\frac d{du}\hbox{|jarray(l,a)|}, \cr
}$$
for $-|lmax|-1 \le l \le |lmax|$ and $0\le a \le |amax|$. The dimensioning
of the arrays is controlled by |lmax1| which should be at least $|lmax|+2$,
in order to provide additional work space.
@ Here is |legendrejn|.
@<Functions...@>=
subroutine legendrejn(u,c,lmax,jarray,djarray,lmax1)
implicit_none
integer j0,j1,j2,j02,j11,j22,j022
parameter (j0=0, j1=1, j2=2, j02=3, j11=4, j22=5, j022=6)
integer l,a,lmax,amax,lmax1
real u,c,jarray,djarray,scale,g
dimension jarray(-lmax1-1:lmax1,0:6),djarray(-lmax1-1:lmax1,0:6)
external legendrejs
amax=2
call legendrejs(u/c,lmax+2,amax,jarray,lmax1)
scale=1
do l=1,lmax+2
scale=scale*u/(2*l+1)
do a=0,amax
jarray(l,a)=scale*jarray(l,a)
end do
end do
scale=1
do l=-1,-lmax-2,-1
scale=scale/u*(2*l+3)
do a=0,amax
jarray(l,a)=scale*jarray(l,a)
end do
end do
do l=-lmax-2,lmax /* Calculate multiple index solutions */
jarray(l,j02) = u/2*jarray(l+1,1)
jarray(l,j11) = u/2*jarray(l+1,0)
jarray(l,j22) = u/2*jarray(l+1,1) + (u/c)^2/2*jarray(l+2,0)
jarray(l,j022) = u^2/8*jarray(l+2,0)
end do
g = sqrt(1 + (u/c)**2) /* Calculate derivatives */
do a=0,6
do l=-lmax-1,lmax
djarray(l,a) = jarray(l-1,a)/g - (l+1)/u*jarray(l,a)
end do
end do
return
end
@* Index.
